function y = powermod_z_m(m, CONSTS)

    r_high = 1e+3;
    
%     inside = quad_gauss(@(r)real(func_power_inside(r, m, CONSTS)), [0, 1], 0.002, 4);
%     outside = quad_gauss(@(r)real(func_power_outside(r, m, CONSTS)), [1, r_high], 0.002, 4);

    inside = quadgk(@(r)real(func_power_inside(r, m, CONSTS)), 0, 1);
    outside = quadgk(@(r)real(func_power_outside(r, m, CONSTS)), 1, r_high);

    
    y = (inside + outside);

end

function y = func_power_inside(r, m, CONSTS)

    a = CONSTS.a;
    c = CONSTS.c;
    a_mn = func_coef_a_mn(m, CONSTS);
    
    [E_rho_inside, E_phi_inside, foo] = func_E_inside(r, m, CONSTS);
    [H_rho_inside, H_phi_inside, foo] = func_H_inside(r, m, CONSTS);
    
    y = (c/4)*a^2*(abs(a_mn))^2*(E_rho_inside.*conj(H_phi_inside) - ...
                     E_phi_inside.*conj(H_rho_inside)).*r;

    y(isnan(y))=0;

end

function y = func_power_outside(r, m, CONSTS)

    a = CONSTS.a;
    c = CONSTS.c;
    a_mn = func_coef_a_mn(m, CONSTS);
    
    [E_rho_outside, E_phi_outside, foo] = func_E_outside(r, m, CONSTS);
    [H_rho_outside, H_phi_outside, foo] = func_H_outside(r, m, CONSTS);
    
    y = (c/4)*a^2*(abs(a_mn))^2*(E_rho_outside.*conj(H_phi_outside) - ...
                     E_phi_outside.*conj(H_rho_outside)).*r;
                 
    y(isnan(y))=0;

end

function y = func_coef_a_mn(m, CONSTS)

    d = CONSTS.d;
    a = CONSTS.a;
    k0 = CONSTS.k0;
    
    N_mn = func_norm(m, CONSTS);
    coef_Im = (1/pi)*func_coef_Im(m, CONSTS);
    if(m==0)
        coef_Im = (1/pi)*func_a0(CONSTS);
    end
    
    [foo, p_mn] = func_coef_of_fields(m, CONSTS);
    [foo, E_phi_m_border, foo]  = func_E_inside(1, m, CONSTS);

    y = ((2*pi*a)/N_mn)*E_phi_m_border*coef_Im*pi*besselj(0, k0*p_mn*d);

end




